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The slip phenomena in thin polymer films confined by either flat or periodically 
corrugated surfaces are investigated by molecular dynamics and continuum simula- 
tions. For atomically flat surfaces and weak wall- fluid interactions, the shear rate 
dependence of the slip length has a distinct local minimum which is followed by 
a rapid increase at higher shear rates. For corrugated surfaces with wavelength 
larger than the radius of gyration of polymer chains, the effective slip length decays 
monotonically with increasing corrugation amplitude. At small amplitudes, this de- 
cay is reproduced accurately by the numerical solution of the Stokes equation with 
constant and rate-dependent local slip length. When the corrugation wavelength is 
comparable to the radius of gyration, the continuum predictions overestimate the 
effective slip length obtained from molecular dynamics simulations. The analysis of 
the conformational properties indicates that polymer chains tend to stretch in the 
direction of shear flow above the crests of the wavy surface. 

PACS numbers: 83.50.Rp, 83.50.Lh, 68.35. Ct, 61.20.Ja, 36.20.Ey 

I. INTRODUCTION 

The dynamics of fluid flow in conflned geometries has gained renewed interest due to 
the recent developments in micro- and nanofluidics l|]. The investigations are motivated 
by important industrial applications including lubrication, coating, and painting processes. 
The flow behavior at the sub-micron scale strongly depends on the boundary conditions at 
the liquid/solid interface. A number of experimental studies on fluid flow past nonwetting 
surfaces have shown that the conditions at the boundary deviate from the no-slip assump- 
tion [2I. The most popular Navier model relates the slip velocity (the relative velocity of the 
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fluid with respect to the adjacent sohd wall) and the shear rate with the proportionality co- 
eflicient, the slip length, which is determined by the linear extrapolation of the fluid velocity 



profile to zero. The magnitude of the slip 



wettability j^, |4j5, 
and shear rate 115 



length 



surface roughness 
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depends on several key parameters, such as 



Id, 
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121], complex fluid structure 



16l . Il7j . However, the experimental determination of the slip length as a 
function of these parameters is hampered by the presence of several factors with competing 
effects on the wall slip, e.g. surface roughness and wettability [3] or surface roughness and 
shear rate 



In recent years, molecular dynamics (MD) simulations have been wide 



the slip flow past atomically smooth, homogeneous surfaces 
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271]. The advantage of the MD approach is that the velocity proflles and shear stresses 



are resolved at the molecular level. The slip length in the shear flow of simple fluids past 
crystalline walls is a function of the wall-fluid density rati o |l9l. 



atoms and fluid molecules 
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261 ] ■ the surface energy [19 



2l| . the relative size of wall 



271], and the interfacial shear 



271 ]. Weak wall- fluid interactions and incommensurable structures of the solid 
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and fluid phases at the interface usually lead to enhancement of slip 
length at low shear rates is about several molecular diameters then it increases with the shear 
.ate. a.d the slope of the .ate dependence ,s gteate. fo. weaUet wa..-fln,d inte.act.ons HQ. 
The rate dependence of the slip length in the flow of polymer melts is more complicated 



because of the additional length and time scales associated with the 



chains at the interface and the shear thinning yiscosity 
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dynamics of polymer 
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In t he p resence of surface roughness 



40|, the fluid flow near the solid boundary is perturbed on the length scales of the 



surface heterogeneities and its description requires deflnitions of the effective slip length and 
the average location of the reference plane. Most commonly, the location of the reference 
plane is deflned as the mean height of the surface asperities, and the shear rate is determined 
by averaging of the fluid flow over the typical length scale of the surface inhomogeneities. 
In general, the surface roughness is expected to reduce the effective slip length for wetting 



liquids 
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351 ]. For sufficiently rough surfaces the no-slip boundary condition can be 



achieved even if the local condition is of zero shear stress 



411 ]. However, in special cases, when 



the fluid is partially dewetted at the nanostructured, or the so-cal 



surfaces, the slip length might be enhanced up to a few microns 42 
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In recent MD studies on shear flow of simple fluids, the behavior of the effective shp length 
was investigated in the Couette cell with either mixed boundary conditions [38] or periodic 
surface roughness 33|. In the flrst study, the lower stationary wall with mixed boundary 
conditions was patterned with a periodic array of stripes representing alternating regions of 
flnite slip and zero shear stress. In the other study [331], the periodically roughened surface 
was modeled by introducing a sinusoidal offset to the position of the wall atoms. At the wavy 
wall, the local slip length is modifled by the presence of curvature and becomes position- 
dependent along the curved boundary Q, 5^. A detailed comparison between continuum 
analysis and MD simulations shows an excellent agreement between the velocity proflles 
and effective slip lengths when the characteristic length scale of substrate inhomogeneities is 
larger than approximately thirty molecular diameters 33|, l38|] . In the case of rough surfaces. 



an additional correction due to variable wall density was incorporated in the analysis [33]. 
The problem of applicability of the results obtained for monoatomic fluids to polymer melts 
is important for modeling polymer flows in confined geometries and design of the hybrid 
continuum-atomistic algorithms. 

In this paper, the MD simulations are carried out to study the dynamic behavior of 
the slip length at the interface between a polymer melt and atomically flat or periodically 
corrugated surfaces. The MD results for flat crystalline walls conflrm previous flndings [30] 
that the slip length goes through a local minimum at low shear rates and then increases 
rapidly at higher shear rates. For periodically corrugated surfaces and wetting conditions, 
the effective slip length decreases gradually with increasing values of the wavenumber. The 
solution of the Stokes equation with either constant or rate-dependent local slip length 
is compared with the MD simulations for corrugated surfaces with wavelengths ranging 
from molecular dimensions to values much larger than the radius of gyration of polymer 
chains. The orientation and the dynamics of linear polymer chains are signiflcantly affected 
by surface roughness when the corrugation wavelengths are comparable with the radius of 
gyration. 

The rest of the paper is organized as follows. The details of molecular dynamics and 
continuum simulations are described in the next section. The dynamic response of the slip 
length and shear viscosity in the cell with atomically flat surfaces is presented in Sec. lIIII The 
results of MD simulations for corrugated walls with large wavelengths and comparison with 
continuum predictions are reported in Sec. lIVI The slip behavior for small wavelengths and 
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the conformational properties of the polymer chains near the rough surfaces are analyzed in 
Sec.|Vl The summary is given in the last section. 

II. THE DETAILS OF THE NUMERICAL SIMULATIONS 

A. Molecular dynamics model 

The computational domain consists of a polymeric fluid confined between two atomistic 
walls. Figure[I] shows the MD simulation setup. Any two fluid monomers within a cut-off 
distance of = 2.5 a interact through the truncated Lennard- Jones (LJ) potential 
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VLj(r) = As 



a\^^ fay 
r ) 



(1) 



where e is the energy scale and a is the length scale of the fluid phase. The LJ potential was 
also employed for the wall- fluid interaction with e^,/ = e and cr^/ = a . The wall atoms do 
not interact with each other, and the wall-fluid parameters are fixed throughout the study. 
In addition to the LJ potential, the neighboring monomers in a polymer chain (A^ = 20 
beads) interact with the finite extensible nonlinear elastic (FENE) potential 



1 



VFENE{r) = ~krl\^ , (2) 

/ L \Tq/- 

where To = 1.5 a and fc = 30ecr~^ fsi)]. The MD simulations were performed at a constant 

density ensemble with p = 0.88(7~^. The total number of fluid monomers is A^/= 67200. 

The motion of the fluid monomers was weakly coupled to an external thermal reser- 
I I 

voir [52] . To avoid a bias in the flow direction, the random force and the friction term were 
added to the equation of motion in the y direction [3] 

where F = 1.0 is a friction constant that regulates the rate of heat flux between the fluid 
and the heat bath, and fi{t) is the random force with zero mean and variance 2mTkBT5{t) 



determined from the fluctuation-dissipation theorem 



531]. Temperature of the Langevin 



thermostat is set to T = l.le/kB, where is the Boltzmann constant. The equations of 



motion were integrated using the Verlet algorithm 5J] with a time step At = 0.005 r, where 
r = ^Jma"^ je is the characteristic LJ time. 
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FIG. 1: (Color online) A snapshot of the fluid monomers (open circles) confined between solid 
walls (closed circles) obtained from the MD simulations. The atoms of the stationary lower wall 
are distributed along the sinusoidal curve with the wavelength A and amplitude a. The flat upper 
wall is moving with a constant velocity U in the x direction. The effective slip length L^s is 
determined by the linear extrapolation of the velocity profile to Ux = 0. 

The dimensions of the system in the xy plane, unless specified otherwise, were set to 
Lx = 66.60 a and Ly = 15.59 cr. The upper wall was composed of two layers of an fee lattice 
with density = 1.94 cr~^, which corresponds to the nearest-neighbor distance of d = 0.9a 
between wall atoms in the (111) plane. The lower wall was constructed of two fee layers 
of atoms distributed along the sinusoidal curve with the wavelength A and amplitude a. 
For the largest wavelength A = 66.60 a, the density of the lower wall pw = 1.94 cr~'^ was kept 
uniform along the sinusoid (by including additional rows of atoms parallel to the y axis) to 
avoid additional analysis of the effective slip length due to variable wall density j^S]- In the 
present study, the corrugation amplitude was varied in the range ^ a/a ^ 12.04. In the 
absence of the imposed corrugation (a = 0) the distance between the inner fee planes is set 
to Lz = 74.15 a in the z direction. Periodic boundary conditions were imposed in the x and 
y directions. 

The initial velocities of the fluid monomers were chosen from the Maxwell-Boltzmann 
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probability distribution at the temperature T =l.le/kB- After an equilibration period of 
about 3 X 10"^ r with stationary walls, the velocity of the upper wall was gradually increased 
in the x direction from zero to its final value during the next 2 x 10'^ r. Then the system 
was equilibrated for an additional period of 6 x 10'^ r to reach steady-state. Averaging time 
varied from 10^ r to 2 x 10^ r for large and small velocities of the upper wall respectively. 
The velocity profiles were averaged within horizontal slices of x Ly x Az, where Az = 
0.2(7. Fluid density profiles near the walls were computed within slices with thickness 
Az = O.Ola [27]. 



B. Continuum method 



A solver based on the finite element method was developed for the two-dimensional 
steady-state and incompressible Navier-Stokes (NS) equation. The NS equation with these 
assumptions is reduced to 

p (u ■ Vu) = - Vp + /iV^u, (4) 

where u is the velocity vector, p is the fluid density, and p and p are the pressure field and 
viscosity of the fluid respectively. 

The incompressibility condition is satisfied by a divergence-free velocity field u. In order 
to avoid the decoupling of the velocity and the pressure fields in the numerical simulation 
of the incompressible flow, the penalty formulation is adopted 551] • This method replaces 
the continuity equation, V ■ u = 0, with a perturbed equation 



Vu 



P 
A' 



(5) 



where A is the penalty parameter. For most practical applications, where computation is 
performed with double-precision 64 bit words, a penalty parameter between 10^ and 10^ is 
sufficient to conserve the accuracy |55|]. In our simulations the incompressibility constraint 
was set to A = 10''. The pressure term in Eq. is then substituted into the NS equation. 
The equation Eq. (jl]) can be rewritten as follows: 



p (u ■ Vu) = AV(V ■ u) + 



(6) 



where the continuity equation is no longer necessary 55] . The NS equation is integrated 



with the Galerkin method using bilinear rectangular isoparametric elements 



55|. 
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Four boundary conditions must be specified for tlie continuum simulation. Periodic 
boundary conditions are used for the inlet and outlet along the x direction. A slip boundary 
condition is apphed at the upper and lower walls. In the local coordinate system (spanned 
by the tangential vector t and the normal vector n), the fluid velocity along the curved 
boundary is calculated as 

ut = Lo[{n-V)ut + Ut/R], (7) 

where Ut is the tangential component of u = utt + UnU, Lq is the slip length at the flat 
liquid/solid interface, the term in the brackets is the local shear rate, and R is the local 
radius of curvature |5(|. The radius of curvature is positive for the concave and negative 
for the convex regions. For a fiat surface, R —>■ oo, the boundary condition given by Eq. ([7]) 
simply becomes the Navier slip law. 

The simulation is started by applying the no-slip boundary condition as the initial guess. 
Once the equations of motion are solved implicitly, the local slip velocities at the lower and 
upper boundaries are updated using Eq. ([7]). In the next step, the equations of motion are 
solved with the updated slip velocities used as a new boundary condition. The iterative 
procedure is repeated until the solution converges to a desired accuracy. The convergence 
rate of the iterative solution remains under control with the under-relaxation value about 
0.001 for the boundary nodes. In all continuum simulations, the grids at the lower boundary 
have an aspect ratio of about one. The computational cost is reduced by increasing the aspect 
ratio of the grids in the bulk region. 

The normalized average error value in the simulation is defined as 



error 

i=l 



where Np is the number of nodes in the system, u" is the velocity at the node i and time 
step n, and u""*"^ is the velocity in the next time step. The typical error in the converged 

;„ 1 j-1 1 n— 9 A j. _ i _ r dut — i _ 

■'local g„ 5 

local slip length is Liocai = (-^o ^ — R 



III. MD RESULTS FOR FLAT WALLS 

The averaged velocity profiles for selected values of the upper wall speed U are presented 
in Fig.[2j The profiles are linear throughout the cell, except for U ^ 6.0 a/r where a slight 
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FIG. 2: (Color online) Averaged velocity profiles in the cell with flat upper and lower walls. The 
solid lines are the best linear fit to the data. The vertical axes indicate the location of the fee 
lattice planes. The velocities of the upper wall are tabulated in the inset. 

curvature appears in the region of about 4 a near the walls. Note that the relative slip 
velocity at the upper and lower walls increases with increasing upper wall speed. The shear 
rate was determined from the linear fit to the velocity profiles across the whole width of 
the channel (see Fig. [2]). The uncertainty in the estimated value of shear rate is due to the 
thermal fluctuations and the slight curvature in the velocity profiles near the walls. The 
typical error bars for the shear rate are about 2 x lO^^r^^ and 6 x 10~^r~^ for small and 
large upper wall speeds respectively (not shown). 

In this study, the shear stress in steady-state flow was computed from the Irving-Kirkwood 
relation j56|. The dynamic response of the fluid viscosity with increasing shear rate is 
presented in Fig.O At higher shear rates, the fluid exhibits shear thinning behavior with 
the slope of about —0.33. Although the power law coefficient is larger than the reported 
values in experimental studies {5'(|, the results are consistent with previous MD simulations 



of polymer melts for similar flow conditions 30|, l58|, |59 |. 



The slip length was calculated by the linear extrapolation of the fluid velocity profile 
to zero with respect to a reference plane, which is located 0.5 a away from the fee lattice 
plane FigureH] shows the rate dependence of the slip length in the same range 

of shear rates as in Fig.O The slip length goes through a shallow minimum at low shear 
rates and then increases rapidly at higher rates. The error bars are larger at low shear rates 
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FIG. 3: (Color online) Viscosity of the polymer melt /x/ero" as a function of shear rate. The 
dashed line with the slope —0.33 is plotted for reference. 

because the thermal fluid velocity = kBT/m is greater than the average flow velocity. 
The nonmonotonic behavior of the slip length in sheared polymer films with atomically flat 
surfaces can be interpreted in terms of the friction coefficient at the liquid/solid interface, 
which undergoes a gradual transition from a nearly constant value to the power law decay 



function of the slip velocity 



30|. The data for the slip length shown in Fig.lHare well 



fitted by the fourth-order polynomial 

Lo{x) /a = 16.8 - 72.0 x 10 x + 44.0 x 10^ - 97.3 x 10^ + 80.5 x 10^ x^ (9) 

where x = 7r is the shear rate. The polynomial fit will be used to specify the boundary 
conditions for the continuum solution described in the next section. 



IV. RESULTS FOR PERIODICALLY CORRUGATED WALLS: LARGE 

WAVELENGTH 

A. MD simulations 

The monomer density profiles were computed in the averaging regions located in the 
grooves and above the peaks of the corrugated lower wall with wavelength A = = 66.60 cr. 
The dimensions of the averaging regions are set to 0.5 cr and Ly = 15.59 a in the x and y 
directions respectively. FigureO shows the monomer density profiles near the upper and 
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FIG. 4: (Color online) Variation of the slip length as a function of shear rate in the cell with flat 
upper and lower walls. The solid line is a fourth-order polynomial fit to the data given by Eq. ([UJ. 

lower walls at equilibrium. The pronounced density oscillations are attributed to the suc- 
cessive layering of the fluid monomers near the walls. These oscillations decay within a few 
molecular diameters from the walls to a uniform proflle characterized by the bulk density of 
p = 0.88 cr"^. Note that the height of the flrst peak in the density proflle inside the groove is 
slightly larger than its value above the crest [see Fig.[5](b)]. The fluid monomers experience 
stronger net surface potential in the groove than above the crest because of the closer spatial 
arrangement of the wall atoms around the location of the flrst density peak in the groove. 
This effect is amplifled when the local radius of curvature at the bottom of the grooves is 
reduced at smaller wavelengths (see below). 

The averaged velocity proflles in the cell with periodically corrugated lower and flat upper 
walls are presented in Fig. El The fluid velocity proflles were averaged within horizontal slices 
of thickness Az = 0.2 cr distributed uniformly from the bottom of the grooves to the upper 
wall. With increasing corrugation amplitude, the slip velocity decreases and the proflles 
acquire a curvature near the lower wall. Note also that the relative slip velocity at the upper 
wall increases with the shear rate. The linear part of the velocity proflles, 30 ^ z/cr ^ 60, was 
used to determine the effective slip length Lcs at the corrugated lower wall. The variation 
of the effective slip length as a function of wavenumber ka = 2Txa/ \ is shown in Fig. [71 The 
slip length decreases monotonically with increasing wavenumber and becomes negative at 
ka > 0.7. The results of comparison between the MD and continuum simulations for three 
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FIG. 5: (Color online) Averaged density profiles near the stationary upper wall (a), above the peak 
and in the groove of the lower wall with amplitude a = 8.16 ex and wavelength A = 66.60 o" (b). 

different cases are presented in the next section. 

B. Comparison between MD and continuum simulations 

In continuum simulations, the length, time and energy scales are normalized by the LJ 
parameters a, r and e respectively. The continuum nondimensional parameters are denoted 
by the (~) sign. The size of the two-dimensional domain is fixed to 66.60 x 73.15 in the x 
and z directions, respectively. The following three cases were examined: the Stokes solution 
with constant slip length in Eq. ([7]), the Stokes solution with shear-rate-dependent slip length 
given by Eq. ([9]), and the Navier-Stokes solution with constant slip length in Eq. ([7]). For all 
cases considered, the fiat upper wall is translated with a constant velocity U = 0.5 in the x 
direction. Similarly to the MD method, the effective slip length is defined as a distance from 
the reference plane at a = to the point where the linearly extrapolated velocity profile 
vanishes. 

In the first case, the finite element method was implemented to solve the Stokes equation 
with boundary conditions at the upper and lower walls specified by Eq. ([71) with Lq = 14.1. 
This value corresponds to the slip length Lq = 14.1 ±0.5 a extracted from the MD simulations 
in the cell with fiat walls and the velocity of the upper wall U = 0.5 a /r. Note that at 
large amplitudes a ^ 8.16, the normal derivative of the tangential velocity dut/dn at the 
bottom of the grooves is negative, while the x component of the slip velocity is positive 
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FIG. 6: Averaged velocity profiles for the indicated values of the corrugation amplitude a. The 
vertical dashed line denotes a reference plane for calculation of the effective slip length at the 
corrugated lower wall. The velocity of the flat upper wall is U= 0.5 a /t. 

everywhere along the corrugated lower wall. The dependence of the effective slip length on 
the corrugation amplitude is shown in Fig. [71 The continuum results agree well with the 
approximate analytical solution 50|] for ka < 0.5 (not shown). For larger amplitudes, ka > 
0.5, where the analytical solution is not valid, our results were tested to be grid independent. 
There is an excellent agreement between slip lengths obtained from the MD and continuum 
simulations for ka < 0.3. With further increasing the amplitude, the slip length obtained 
from the continuum solution overestimates its MD value. The results presented in Fig. [7] 
are consistent with the analysis performed earlier for simple fluids js^, although a better 
agreement between MD and continuum solutions was expected at ka > 0.5 because of the 
larger system size considered in the present study. 

As discussed in the previous section, the slip length for atomically flat walls is rate- 
dependent even at low shear rates (see Fig.H]). In the second case, we include the effect of 
shear rate in the analysis of the effective slip length at the corrugated lower wall and flat 
upper wall. The Stokes equation is solved with boundary conditions given by Eq. ([7]), where 
the slip length Eq. ([9]) is a function of the local shear rate at the curved and flat boundaries. 
The results obtained from the Stokes solution with constant and rate-dependent slip lengths 
are almost indistinguishable (see Fig. [7]). This behavior can be attributed to a small variation 
of the intrinsic slip length Eq. ([9]) at low shear rates. For example, at the largest amphtude. 
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FIG. 7: The effective slip length as a function of wavenumber ka obtained from the MD simulations 
(♦), the solution of the Stokes equation with rate- independent slip length Lq (o) and with Lq{'^t) 
given by Eq. ([9]) (x), the solution of the Navier-Stokes equation with Lq (a). 

a = 12.04, the local shear rate at the corrugated wall is position-dependent and bounded by 
\dut/dn + Ut/R\ ^ 0.0035. In this range of shear rates, the normalized value of the slip 
length in Eq. ([9]) varies between 14.7 ^ Lq ^ 16.6. It is expected, however, that the effect 
of shear rate will be noticeable at larger values of the top wall speed U . 

In the third case, the Navier-Stokes equation is solved with a constant slip length Lq = 14.1 
in Eq. ([7]) at the flat upper and corrugated lower walls. The upper estimate of the Reynolds 
numbe. based on the BuM density p = 0.88, v.co.ity A = 20.0, a.d the flu.d velo».ence 
across the channel is Re ~ 1.3. It was previously shown by Tuck and Kouzoubov [61] that at 
small ka and Lo = the magnitude of the apparent slip velocity at the mean surface increases 
due to finite Reynolds number effects for i?e > 30. In our study, the difference between the 
slip lengths extracted from the Stokes and Navier-Stokes solutions is within the error bars 
(see Fig. [7]). These results confirm that the slip length is not affected by the inertia term in 
the Navier-Stokes equation for Re < 1.3. To check how sensitive the boundary conditions 
are to higher Reynolds number flows, we have also repeated the continuum simulations for 
larger velocity of the upper wall U= 50, which corresponds to i?e ~ 130. For the largest 
corrugation amplitude a = 12.04, the backfiow appears inside the groove and the effective 
slip length becomes smaller than its value for U = 0.5 by about 0.7 (not shown). 
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FIG. 8: Averaged fluid density profiles near Hat upper wall (a) and corrugated lower wall with 
amplitude a = lAa and wavelength A = 7.5(7 (b). The velocity of the upper wall is U = 0.5 a /t. 

V. RESULTS FOR PERIODICALLY CORRUGATED WALLS: SMALL 

WAVELENGTHS 



A. Comparison between MD and continuum simulations 

The MD simulations described in this section were performed at corrugation wavelengths 
(A/cr = 3.75, 7.5 and 22.5) comparable with the size of a polymer coil. Periodic surface rough- 
ness of the lower wall was created by displacing the fee wall atoms by Az = a sin (2 ttx/X) in 



the z direction [33|| • In order to reduce the computational time the system size was restricted 
to Nf = 8580 fluid monomers and = 22.5 a, Ly = 12.5 a and Lz = 35.6 a. All other system 
parameters were kept the same as in the previous section. 

The representative density proflles near the upper and lower walls are shown in Fig. [8] for 
the wavelength A = 7.5 a. The height of the flrst peak in the density proflle is larger in the 
grooves than near the flat wall or above the crests of the corrugated surface. The effective 
slip length as a function of wavenumber ka is plotted in Fig.O For all wavelengths, the slip 
length decreases monotonically with increasing values of ka. At the smallest wavelength 
A = 3.75(7, the slip length rapidly decays to zero at A;a~0.4 and weakly depends on the 
corrugation amplitude at larger ka. Inspection of the local velocity proflles for A = 3.75 a 



and ka > 1.0 indicates that the flow is stagnant inside the grooves 60 1 



In the continuum analysis, the Stokes equation with a constant slip length Lq = 14.1 in 
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FIG. 9: The effective slip lengtli as a function of wavenumber ka for the indicated values of 
wavelength A. Continuum results are denoted by the dashed lines and open symbols, while the MD 
results are shown by straight lines and filled symbols. The inset shows the same data for ka ^ 0.25. 

Eq. (171) is solved for the three wavelengths. The comparison between the MD results and 
the solution of the Stokes equation is presented in Fig. [91 The error bars are larger for the 
smallest wavelength because of the fine grid resolution required near the lower boundary at 
ka < 0.2 (see inset in Fig.[9|). The results shown in Fig. [9l confirm previous findings for simple 
fluids 331] that the slip length obtained from the Stokes flow solution overestimates its MD 
value and the agreement between the two solutions becomes worse at smaller wavelengths. 
It is interesting to note that the curves for different wavelengths intersect each other at 
ka ~ 0.63 in the MD model and at ka ~ 1.02 in the continuum analysis. The same trend was 
also observed in the previous study on slip flow of simple fluids past periodically corrugated 
surfaces [331]. 



B. The polymer chain configuration and dynamics near rough surfaces 

In this section, the properties of polymer chains are examined in the bulk and near the 
corrugated boundary with wavelengths A = 3.5o" and A = 7.5cr. The radius of gyration Rg 
was computed as 



TV 

32 ' 



Rt = j^^2^R^'Rcm)^ (10) 

1=1 
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A = 7.5cj 


a = 0.2cr 


a = 0.6 a 


a = 1.4cr 




Rgx 


^gy 


Rgz 


Rgx 


^gy 


Rgz 


Rgx 


^gy 


Rgz 


Bulk Equilibrium 


1.18 


1.18 


1.18 


1.18 


1.18 


1.18 


1.18 


1.18 


1.18 


Shear flow 


1.70 


1.11 


1.04 


1.76 


1.10 


1.02 


1.79 


1.10 


1.01 


Upper Equilibrium 


1.34 


1.37 


0.66 


1.35 


1.36 


0.66 


1.35 


1.35 


0.66 


wall Shear flow 


1.63 


1.29 


0.63 


1.64 


1.29 


0.63 


1.67 


1.29 


0.63 


Peak Equilibrium 


1.35 


1.36 


0.69 


1.40 


1.31 


0.76 


1.46 


1.27 


0.84 


Shear flow 


1.65 


1.28 


0.67 


1.85 


1.21 


0.75 


2.17 


1.10 


0.93 


Groove Equilibrium 


1.31 


1.37 


0.63 


1.10 


1.47 


0.61 


0.77 


1.69 


0.61 


Shear flow 


1.50 


1.34 


0.62 


1.19 


1.44 


0.60 


0.75 


1.89 


0.59 



TABLE I: Averaged x, y, and z components of the radius of gyration at equilibrium and in the 
shear flow. The Rg values are reported in the bulk, near the flat upper wall, above the peaks, and 
inside the grooves. The wavelength of the lower wall is A = 7.5 a. The size of the averaging region 
inside the grooves and above the peaks is cr x 12.5 o" x 1.5 cr. The estimate of the error bars is 
±0.03(7. 

where Ri is the three-dimensional position vector of a monomer, = 20 is the number of 
monomers in the chain, and Rem is center of mass vector defined as 

1 ^ 

Rcm = j;^^Ri- (11) 

4 = 1 

The chain statistics were collected in four different regions at equilibrium {U = 0) and in the 
shear flow induced by the upper wall moving with velocity U = 0.5 a /t in the x direction. 
Averaging regions were located above the peaks, in the grooves, near the flat upper wall and 
in the bulk (see Fig. [10] for an example). The dimensions of the averaging regions above the 
peaks and in the grooves of the lower wall are a x 12.5 cr x 1.5 cr, and near the upper wall 
and in the bulk are 22.5 cr x 12.5 a x 1.5 cr. Three components of the radius of gyration were 
computed for polymer chains with the center of mass inside the averaging regions. 

In the bulk region, the components of the radius of gyration remain the same for both 
wavelengths, indicating that the chain orientation is isotropic at equilibrium and is not af- 
fected by the confining walls. In the steady-state flow, the effective slip length is suppressed 
by the surface roughness and, therefore, the shear rate in the bulk increases with the corru- 
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A = 3.75 a 


a = 0.07 a 


a = 0.2o- 


a = l.Oa 




Rgx 




Rgz 


Rgx 


^gy 


Rgz 


Rgx 


^gy 


Rgz 


Bulk Equilibrium 


1.18 


1.18 


1.18 


1.18 


1.18 


1.18 


1.18 


1.18 


1.18 


Shear flow 


1.70 


1.12 


1.03 


1.76 


1.10 


1.02 


1.78 


1.10 


1.01 


Upper Equilibrium 


1.33 


1.36 


0.66 


1.36 


1.34 


0.66 


1.37 


1.34 


0.65 


wall Shear flow 


1.60 


1.32 


0.63 


1.64 


1.29 


0.63 


1.65 


1.27 


0.63 


Peak Equilibrium 


1.38 


1.33 


0.67 


1.42 


1.30 


0.71 


1.41 


1.25 


0.92 


Shear flow 


1.61 


1.27 


0.65 


1.66 


1.26 


0.69 


1.78 


1.17 


0.84 


Groove Equilibrium 


1.33 


1.36 


0.64 


1.25 


1.42 


0.61 


0.48 


2.64 


0.55 


Shear flow 


1.60 


1.29 


0.61 


1.60 


1.33 


0.59 


0.47 


2.70 


0.54 



TABLE II: Averaged x, y, and z components of the radius of gyration at equilibrium and in the 
shear flow. The Rg values are reported in the bulk, near the flat upper wall, above the peaks, and 
inside the grooves. The wavelength of the lower wall is A = 3.75 a. The dimensions of the averaging 
region are the same as in TableHl 

gation amplitude. This explains why the x component of the radius of gyration Rg^ increases 
slightly at larger amplitudes (see Tables!!] and [IT|) . Near the upper wall, the polymer chains 
become flattened parallel to the surface and slightly stretched in the presence of shear flow. 
These results are consistent with the previous MD simulations of polymer melts confined 



between atomically fiat walls 



28 



62 



63|. 



In the case of a rough surface with the wavelength A = 7.5 a, a polymer chain can be 
accommodated inside a groove (see Fig. [10] for an example). With increasing corrugation 
amplitude, the polymer chains inside the grooves elongate along the y direction and contract 
in the x direction (see Table[T]). The tendency of the trapped molecules to orient parallel to 



the grooves was observed previously in MD simulations of hexadecane [32| . In the presence 
of shear flow, polymer chains are highly stretched in the x direction above the crests of the 
wavy wall. A snapshot of the unfolded chains during migration between neighboring valleys 
is shown in Fig.[TOl The flow conditions in Fig. [10] correspond to a negative effective slip 
length Leff ~ —2 a. 

For the smallest corrugation wavelength A = 3.75 a, polymer chains cannot easily fit in 
the grooves unless highly stretched. Therefore, the y component of the radius of gyration 
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l.So 




FIG. 10: (Color online) A snapshot of four polymer chains near the lower corrugated wall for 
wavelength A = 7.5 cr and amplitude a = 1.4(7. The velocity of the upper wall is U = 0.5 u/r. 

is relatively large when the center of mass is located in the deep grooves (see TablellTll . 
Figure[TT] shows a snapshot of several polymer chains in contact with the lower corrugated 
wall. The chain segments are oriented parallel to the grooves and stretched above the crests 
of the surface corrugation. Visual inspection of the consecutive snapshots reveals that the 
chains near the corrugated wall move, on average, in the direction of shear flow; however, 
their tails can be trapped for a long time because of the strong net surface potential inside 
the grooves. For large wavenumbers ka > 0.5, the magnitude of the negative effective 
slip length is approximately equal to the sum of the corrugation amplitude and Rgz of the 
polymer chains above the crests of the wavy wall. 



VI. CONCLUSIONS 



In this paper the effects of the shear rate and surface roughness on slip flow of a poly- 
mer melt was studied using molecular dynamics and continuum simulations. The linear 
part of the velocity proflles in the steady-state flow was used to calculate the effective slip 
length and shear rate. For atomically flat walls, the slip length passes through a shallow 
minimum at low shear rates and then increases rapidly at higher shear rates. In the case 
of periodic surface heterogeneities with the wavelength larger than the radius of gyration, 
the effective slip length decays monotonically with increasing the corrugation amplitude. 
For small wavenumbers, the effective slip length obtained from the solution of the Stokes 
equation with constant and shear-dependent local slip length is in a good agreement with 
its values coraputed from the MD simulations, in accordance with the previous analysis for 
simple fluids 33|. At low Reynolds numbers, the inertial effects on slip boundary conditions 
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FIG. 11: (Color online) A snapshot of five polymer chains near the lower corrugated wall for the 
wavelength A = 3.75 cr and amplitude a = 1.0a. The figure shows the side view (top) and the top 
view (bottom). The velocity of the upper wall is U = O.Scr/r. 

are negligible. When the corrugation wavelengths are comparable to the radius of gyration, 
polymer chains stretch in the direction of shear flow above the crests of the surface corru- 
gation, while the chains located in the grooves elongate perpendicular to the flow. In this 
regime the continuum approach fails to describe accurately the rapid decay of the effective 
slip length with increasing wavenumber. 
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